{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import signal\n",
    "import numpy as np\n",
    "from numpy import pi, diff, unwrap, angle\n",
    "\n",
    "b = [ 0.00782021,  0.,         -0.01564042,  0.,          0.00782021]\n",
    "a = [ 1.,         -3.66102926,  5.09866292, -3.20227714,  0.7660066 ]\n",
    "w, h = signal.freqs(b, a)\n",
    "group_delay1 = -diff(unwrap(angle(h))) / diff(w)\n",
    "group_delay2 = -np.gradient(angle(h), w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import butter, filtfilt\n",
    "order = 2\n",
    "w1 = 0.1\n",
    "b1, a1 = butter(order, w1, btype='low')\n",
    "print('b1:', b1)\n",
    "print('a1:', a1)\n",
    "w2 = [0.04, 0.1]\n",
    "b2, a2 = butter(order, w2, btype='band')\n",
    "print('b2:', b2)\n",
    "print('a2:',a2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import signal\n",
    "#b, a = signal.iirdesign(0.1, 0.3, 5, 50, ftype='cheby1')\n",
    "w1, gd1 = signal.group_delay((b1, a1))\n",
    "w2, gd2 = signal.group_delay((b2, a2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.figure(2, dpi=800)\n",
    "plt.title('Digital filter group delay')\n",
    "plt.plot(w1, gd1, linewidth=0.5, label='group delay 1')\n",
    "plt.plot(w2, gd2, linewidth=0.5, label='group delay 2')\n",
    "plt.ylabel('Group delay [samples]')\n",
    "plt.xlabel('Frequency [rad/sample]')\n",
    "plt.grid(True)\n",
    "plt.legend(loc='upper right') \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, matplotlib.pyplot as plt\n",
    "from scipy import signal\n",
    "\n",
    "fc = 0.2\n",
    "b, a = signal.butter(15, fc)\n",
    "n = np.arange(500)\n",
    "c = 0.0005\n",
    "x = np.exp(-c*(n-250) ** 2)\n",
    "\n",
    "H, w = signal.freqz(b, a, 4096)\n",
    "W, gd = signal.group_delay((b, a), 4096)\n",
    "\n",
    "w0 = .92 * np.pi * fc  # carrier frequency\n",
    "y = x * np.cos(w0*n)   # modulated signal\n",
    "z = signal.lfilter(b,a,y)\n",
    "\n",
    "I = np.argmin([abs(ww-w0) for ww in W])\n",
    "tau = int(gd[I])       # tau\n",
    "\n",
    "plt.plot(W, gd)\n",
    "plt.show()\n",
    "\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(n,y)\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(n,z)\n",
    "plt.plot(n[:-tau],z[tau:], '--')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, matplotlib.pyplot as plt\n",
    "from scipy import signal\n",
    "\n",
    "fc = 0.02\n",
    "fc1, fc2 = 0.019, 0.021\n",
    "\n",
    "b, a = signal.butter(2, [fc1, fc2], btype='bandpass')\n",
    "n = np.arange(5000)\n",
    "x = (abs(n-2000) < 100)\n",
    "\n",
    "H, w = signal.freqz(b, a, 50000)\n",
    "W, gd = signal.group_delay((b, a), 50000)\n",
    "\n",
    "w0 = 1.0 * np.pi * fc  \n",
    "y = x * np.cos(w0*n)  \n",
    "z = signal.lfilter(b,a,y)\n",
    "\n",
    "I = np.abs(W-w0).argmin()\n",
    "tau = int(gd[I]) \n",
    "if tau == 0:            # this happens when scipy has problem to compute group delay and puts 0\n",
    "    tau = int(max(gd))  # then we use the max of group delay, so the shifting \n",
    "                        # should sometimes be *too much*, but never *not enough*\n",
    "\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(n,y)\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(n,z)\n",
    "plt.plot(n[:-tau],z[tau:], '--')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
